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Abstract 

A complete description of the multitudinous ways in which quantum particles can be entangled 
requires the use of high-dimensional abstract mathematical spaces. We report here a particu- 
larly interesting feature of the nine-dimensional convex set-endowed with Hilbert- Schmidt (Eu- 






qh| clidean/flat) measure-composed of two-rebit (4 x 4) density matrices (p). To each p is assigned the 

product of its (nonnegative) determinant \p\ and the determinant of its partial transpose |/5 |- 
negative values of which, by the results of Peres and Horodecki, signify the entanglement of p. 
Integrating this product, |/9||/>^-^| = \pp^'^\, over the nine-dimensional space, using the indicated 
(HS) measure, we obtain the result zero. The two determinants, thus, form a pair of multivariate 
orthogonal polynomials with respect to HS measure. However, orthogonality does not hold, we 
find, if the symmetry of the nine-dimensional two-rebit scenario is broken slightly, nor with the use 
of non-flat measures, such as the prominent Bures (minimal monotone) measure-nor in the full HS 
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j^ ■ extension to the 15-dimensional convex set of two-qubit density matrices. We discuss relations- 

involving the HS moments of \p \-to the long-standing problem of determining the probability 
that a generic pair of rebits/qubits is separable. 
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I. INTRODUCTION AND STATEMENT OF MAIN RESULTS 

The simplest form of finite-dimensional system capable of exhibiting the fundamental, 
holistic property of entanglement-" that feature of quantum formalism which makes it im- 



fl- 



possible to simulate quantum correlations within any classical formalism" [l[-is that com- 
posed of a pair of two-level systems (quantum bits or "qubits") [2:]. The joint state of two 
qubits is describable (2 x 2 = 4) by a 4 x 4 density matrix p, that is a Hermitian matrix 
(its transpose equalling its complex conjugate) having its four nonnegative diagonal entries 
(probabilities)-as well as its four nonnegative eigenvalues-summing to one |3|. The entirety 
of such 4x4 density matrices with their entries restricted to real values-which will be our 
principal subject of analysis here-forms a nine- dimensional convex set. We note that 2x2 
density matrices with real entries have been termed "rebits" 4j. 

Endowing the nine-dimensional set of pairs of rebits with Hilbert-Schmidt (HS) (Eu- 
clidean/fiat) measure [2|, sec. 14.3] [y, |6|, we arrive at the following trio of computational 
results [7|. The three are: (1) the average value or mean (first moment) of the (necessarily 
nonnegative) determinant of p, denoted \p\ (= det(p)), is 2^ = (2*^ ■ 11 ■ 13)"^; (2) the mean 
of the determinant \p^'^\ of the partial transpose (PT) of p-negative values of which, by the 
celebrated results of Peres [8!] and the Horodecki family [9J , are fully equivalent (for qubit- 
qubit systems-as well as qubit-qutrit systems, representable by 6 x 6 density matrices) to 
the entanglement/nonseparability of p-is —gig = — (2-3-11- 13)"^; and (3) the mean of the 
product of these two determinants, that is IpHp"^^! (= \pp^^\ by the Cauchy-Binet [Gram's] 
Theorem), is zero. (The partial transpose-identifiable with the antiunitary operation of 
time- reversal [lO| (cf. ll|, |l2|)-of a 4 x 4 matrix can be obtained simply by transposing in 



place its four 2x2 blocks.) 

So, by assigning to each generic real 4x4 density matrix p, a value equal to this product 
of determinants, we find that the total (but yet undetermined, necessarily negative) value 
allotted to the entangled states exactly cancels the (necessarily nonnegative) value allotted 
to the disentangled/separable states. This result appears to be special to the HS measure, 
as the analogous trio of results using the well-known Bures (minimal monotone) measure 
2, 



. 14.4] pj is gj^ = 2-13 and, numerically, -0.0030959720 and -1.124478 • 10^^ 
7|. It is a natural hypothesis-but one apparently presenting substantial computational 
challenges to test-that this zero-mean (HS) result extends to the standard 15-dimensional 



convex set of 4 x 4 density matrices with arbitrary com])/ea; off-diagonal entries [13|, sec. 2]. 
(From a "physical point of view one should in general consider the entire set of complex 
density matrices" |5|. We can not, however, unqualifiedly, straightforwardly extend the 
interpretation of the zero-product-mean hypothesis to bipartite quantum systems describable 
by density matrices of dimensions greater than four. Then, p^^ can have more than one 
negative eigenvalue, with an even number of such eigenvalues yielding a positive value for 
the determinant \p^'^\ (cf. [3]). Thus, a nonnegative value of Ip'^^l would not necessarily 
indicate that p is separable.) 

The zero-product-mean HS result does not hold, in general, however, for arbitrary two- 
qubit scenarios. This can be seen by slightly modifying in each of three ways, our basic 
nine-dimensional two-rebit example-t hereby diminishing its inherent symmetry. Firstly, if 
we set one off-diagonal pair of entries (p34 = P43, say) of p to zero, giving us an 8-dimensional 
convex set of density matrices, we obtain a mean of the probability distribution over |p| of 
^ = (2^ ■ 33 ■ ll)-i ^ 0.000210438, for |p^^|, a mean of -^ = -13(2^ ■ 3^ ■ ll)-i ^ 
—0.00136785, and for the product of the two determinants, a nonzero mean of 



2196480 

(2^'' ■ 3 ■ 5 ■ 11 ■ 13)^^ ~ 4.55274 ■ 10^^. Secondly, if we increase the dimensionality from nine 
to ten by letting P34 and P43 be arbitrary complex conjugates of each other, then, the trio 
of means, numerically, is 0.000412154, -0.00082468, and (small, but apparently nonzero) 
3.6035 ■ 10~^. Thirdly, for the eight-dimensional set of 4 x 4 density matrices with real 
entries consisting of those minimally degenerate matrices having one eigenvalue zero-which 



constitute the boundary of the basic nine-dimensional set [15|-the mean of the product of 
the three nonzero eigenvalues is ^ ~ 0.0151515, the mean of Ip''^"^! is —tS^ ~ —0.00210438, 
while the mean of the product of \p^'^\ and the three nonzero eigenvalues is also nonzero, 
that is, — 47I20 ~ —0.0000210438. (Of course, in this last minimally degenerate analysis, the 
determinant |p|, equalhng the product of its /ottr eigenvalues, is zero. So, we replace it, for 
our purposes, by the product of the three generically nonzero eigenvalues.) We, additionally, 
have found that no pair of principal (4 x 4, 3 x 3, 2 x 2) minors, other than \p\ and |p^"^|-one 
minor of p and the other of p^^-are orthogonal over the nine-dimensional set with respect 
to HS measure. 

Since both \p\ and \p^'^\ in the generic (two-rebit) nine-dimensional case are polynomials 
in nine variables, it is suggestive, at least, to consider these two determinants as a pair of 



multivariate orthogonal polynomials (MOPS) 16l-ll8| with respect to the HS measure. Grif- 



fiths and Spano have reviewed multivariate orthogonal polynomials with respect to weight 
measures_given by Dirichlet distributions over simplices (such as that to be given below 
by fl26p [6] over the 3-dimensional simplex formed by the diagonal entries of 4 x 4 density 
matrices [18]). 

Let us also note, in our basic nine-dimensional scenario, that the expected value with 
respect to the Hilbert-Schmidt measure of the determinant of the commutator (appearing 
in the von Neumann equation for the time evolution of p), [p, p^^] = pp^^ — p^'^p, is zero 
too. (Jarlskog found the the determinant of the commutator of mass matrices vanishes "if 

n 

and only if there is no CP nonconservation" [19[.) The analogous result was nonzero, that is 
27675648 ~ 2.8545 • 10^^, wlicu, oucc again, we broke the symmetry by setting ^34 = P43 = 0. 
The range of possible values of \p\ is [0,2~^], that of \p^'^\ is [— 2~^,2^^], while that of 
the product |/9||/9^-^|, in the nine-dimensional scenario, we find, is [— 2~^^ ■ 3^'^, 2^^^]. This 
last upper limit is reached, clearly, by the fully mixed (classical) state-the diagonal density 
matrix with diagonal entries (probabilities) all equal to |. An instance of a density matrix 
attaining the last lower limit is 



6V2 3 12 ( 1 V 3) g^ 

^(-1-73) 



6V2 12 V V / 3 6^2 



(1) 



6V2 6^2 

Now, the determinant \p\ is ^ (2^3-3) ^ 0.000805732 and \p^^\ = ^ (-3-273) ^ 
-0.0112224 (their product being the indicated lower limit -j^^ ~ -9.04225-10"^). Both p 
and p^^ have three identical eigenvalues (^^ (3 — y/S) ~ 0.105662 for p, and ^^ (3 + v^) ~ 
0.394338 for p^^. The isolated eigenvalue for p is i (l + \/3) ^ 0.683013, and for p^^, 
i (1 — y/E) ~ —0.183013. The purity (index of coincidence |2|, p. 56]), Tr(p^), of ([1]) equals 
|, so its inverse (the participation ratio |2|, |20|) is 2. 

A. Failure of Hilbert-Schmidt orthogonality in generic IS-dimensional two-qubit 



case 



With the considerable computational (Mathematica) assistance of Michael Trott, we in- 
vestigated the orthogonality hypothesis in the more demanding 15-dimensional generic two- 
qubit scenario, utilizing the Euler-angle parameterization of Tilnia, Byrd and Sudarshan 



57l |. together with the appropriate formulas of Zyczkowski and Sommers [5|, eq. (3.11)]. For 



U eq- 



the expected value of |p||p-^"^|, we obtained the nonzero value of — 457I264 ~ —2.18519 • 10^''', 
and for the mean of Ip'^^l, the more negative value, —3^ ~ —0.00180599. (Normalizing 
the two quantitates \p\ and \p^'^\ to have unit length [the latter numerically], we find that 
their inner product is approximately -0.07. This translates into an angle of 94 degrees, close 
to orthogonality.) are, since it appears to be even more computationally challenging to nor- 
malize \p^'^\ so that its average HS length is unity.) (However, we are not presently able to 
evaluate how close to orthogonality the determinants \p\ and \p^'^\ are, since it appears to be 
even more computationally challenging to normalize |p^^| so that its average HS length is 
unity.) The mean of IpI is known (a, eq. (3.11)] to be gg^ = (22.3- 17- 19)-^ ^0.000257998. 

B. Separability probabilities 



It has been a relatively long-standing problem of "philosophical, practica 



significance [20| (in the recently-burgeoning field of quantum information 2l|) to determine 



and physical" 



the probability (with respect to a number of measures, such as the HS and Bures) that a 



pair of qubits is separable [20 



22 



251 



31| . In this regard, it would be of interest to, in some 



manner, combine (to achieve a fuller understanding of the "geometry of quantum states" 
[2|) the trio of HS (first) moment (mean) results given above with two theorems (apparently 
the only ones yet developed) pertaining to entanglement in terms of the Hilbert-Schmidt 



measure (cf. 



3l|). 



1. Existing theorems 

One of these theorems states that the probability that a generic bipartite state (of ar- 
bitrary dimension) has a positive partial transpose (PPT) is twice the probability that a 
generic boundary (minimally degenerate [one eigenvalue zero]) state has a PPT. (For the 
proof, it was established that the convex set of mixed PPT states is "pyramid-decomposable 



and hence is a body of constant height" 15|.) The other (but now dimension-specific) 
theorem-derived by enforcing the nonnegativity of pairs of principal 3x3 minors of p^^-is 
that the probability that a generic real 4x4 density matrix has a PPT (or, equivalently, 
by the Peres- Horodecki results, is separable) is no greater than ^^ ^ 0.537619, nor no less 



(using the concept of absolute separability [23| and the important Verstraete-Audenaert-de 
Moor bound [2^) than 6928-2205^ ^ 0.0348338 |25J. (An absolutely separable state is one 
that can not be entangled by any unitary transformation.) From these two theorems to- 
gether, one easily obtains the corollary that the HS separability probability of a generic 
two-rebit minimally degenerate state is no greater than ^^ ~ 0.26881. 



C. Higher-order moments 

In addition to the trio of HS means [first moments) given above, we have been able to 
compute exactly several higher-order moments, as well, for the nine-dimensional scenario. 
The corresponding trio of second (raw/non-central) moments is 2439344 = (2^°-ll-13-17)^^ ~ 
4.01712 ■ 10-', ttJLt = nin,f,... ^ 0.0000108462, and ^ - ^ 



2489344 210-11-13-17 w.wwww j-w"t:w^, ain^ 5696343244800 218-32-52-13-17-19-23 

1.2288585- 10^^^. From our several results, we can deduce that the generic two-rebit Hilbert- 
Schmidt correlation between \p\ and \p^'^\ is positive, that is ,^'^'^ ~ 0.360291. 

Further, making use of the generalized normalization constants |5|, eq. (4.3)], we obtain 
that the m-th moment of the probability distribution over \p\ is 

, _ ^^^ 4^-^-r(2m + 2)r(2m + 4) 

(.m/real-y^^ r(4m + 10) ' ^^ 

Let us also note that the counterpart of this result for the 15-dimensional generic complex 
two-qubit states is 

T{m + l)r(m + 2)r(m + 3)r(m + 4) 
r(4(m + 4)) 



Cm/comple. = 108972864000^ ' \^^^!^ ) ^^^ ' ' L (3) 



1. Application of Chebyshev inequality 

The skewness (71) of the Hilbert-Schmidt probability distribution over Ip"^^] is negative 
(as well as all odd moments [m = 1,3,5,7,9] so far computed), that is, -3.13228-so, the 
left tail of the distribution is more pronounced than the right tail-while its kurtosis (72), 
a measure of "peakedness" is quite high, 17.6316. (Higher kurtosis indicates that more 
of the variance is the result of infrequent extreme deviations than frequent modestly sized 
deviations.) From the first two moments, we obtain the variance 

2 ^ 30397 ^ 9 487838 . ig-e. (4) 

3203785728 ^ ^ 



Application of the standard-form one-sided Chebyshev inequality 26|] (we perform a linear 
transformation, so that negative values of \p^'^\ are mapped to [0,1]), then, yields an upper 
bound on the Hilbert-Schmidt separability probability of the two-rebit density matrices 
of Irfll ~ 0.874759. This is a substantially weaker upper bound, however, than that of 



34749 

ii§ ^ 0.537619 estabhshed in 



25| , by enforcing the nonnegativity of pairs of 3 x 3 principal 



minors of the partial transpose, and even weaker than . " ^ ^ 0.76854, obtained by requiring 



the nonnegativity of all six 2x2 principal minors of the partial transpose 25] 



D. Fits of probability distributions to computed moments 

1. Beta distribution fit to first two moments 

Further, we linearly mapped \p^'^\ E [—jq, ^] to y, so that y G [0, 1], and transformed its 
exact first nine moments accordingly. Then we exactly fit the first and second such moments 
to a basic (two-parameter) beta distribution, giving us a probability distribution of the form 

(Fig.ia, 



v"-''i]-f-'\ „ . i5HH5! « 29.3588. b = 55HH 

B{a,b) ' 516749 ' 2066996 



P\pHy) = - ^T^ ; «=^^T7FiSr^ 29.3588, 6 =55^7^ ^2.42768^ (5) 



where B{a, b) is the beta function. The ratios of the next six transformed moments 
(?7i = 3, . . . , 8) to the corresponding moments of this beta distribution all rather remarkably 
lie between 0.99 and 1 (with the ratio of the ninth moments being 0.986) (Fig. [2]). Integrat- 
ing the distribution ([5]) over the interval y E [if)l]; we obtain an associated separability 
probability estimate of 0.4183149. 

We did similarly linearly transform the product of \p\ and \p^'^\ to lie in the unit interval 
[0,1], and exactly fit to the so-transformed first two moments, a beta distribution of the form 

(cf. (ED) 

y(.a-i)/i_yYb-i) 2392921 , 21536289 

R II pr, y = ^ \ ;' ; a = ^41.4057,6= ^68.8722. 6 

'^"'' '^^^ B{a,b) ' 57792 ' 308224 ^ ' 

The associated probability (over the separability domain [^,1]) was computed as 
0.49331935. Using numerical methods, we found that the (transformed) third moment of the 
distribution of the product |p||p^"^| was approximately 99.84% as large as the third moment 
of (EI), but the fourth moment was only 62.25% as large. So, this fit to a beta distribution 



certainly appears to be inferior to the earlier one 
first two moments of |p'^'^|. 



(Figs. [T] and [2]) based on fitting the 



2. Use of two general moment-fitting procedures 

One can use the calculated exact nine moments of \p^'^\ in certain formulas for cumu- 
lative distribution functions, in order to approximate the desired, specific (separability) 



probability that \p \ is greater than zero [27|, l28|, eq. (2)] (cf. J30|). In Fig. Owe display 



together one such set of nine (greater) separability probability estimates based on the non- 
parametric (stable approximant) reconstruction approach of Mnatsakanov [27|, along with 
the analogous first nine (lesser) estimates using the orthogonal-polynomial-based method- 
ology of Provost and Ha 28|, eq. (3.5)]. (In the latter approach, we used the well-fitting 
beta distribution ([S]) as a "baseline" density that is adjusted by associated modified Jacobi 
orthogonal polynomials.) It is clear that additional (exact) moments (m > 9) are certainly 
needed to more satisfactorily sharpen the generic two-rebit separability probability estimates 
in Fig. [3l Then, one might be able to obtain convincing evidence-in the absence of the 
desired, but highly challenging exact symbolic computation-for the true (hypothetically ex- 



act 



22 



25 



3l| ) value of the separability probability. (Perhaps, somewhat relatedly, Gurvits 



has shown "that the weak membership problem for the convex set of separable normalized 
bipartite density matrices is NP-hard" 32|.) 



II. METHODOLOGIES EMPLOYED 



A. Density Matrix Parameterizations 



Our analysis proceeds in the framework of the Bloore (or correlation coefficient) param- 



eterization 



29 



33 



34| of the 4x4 density matrices (p) which allows us (in the generic 



two-rebit case of immediate interest here) to work primarily in seven dimensions, rather 
than the nine naively expected. Also, in our computations, we still further reparametrize 
three (^13, ^14, Z24,) of the six correlations 



Pij 



y/PiiPjj 



1 < « < i < 4, 



Zij G 



-1,1] 



(7) 



in terms of partial correlations |3J], allowing certain requisite integrations to be per: 



"ormed 



simply over six-dimensional hypercubes, rather than more complicated domains 25|. (Al- 



ternatively, and reasonably computationally competitively, one may utilize the cylindrical 
algebraic decomposition 35| to define the integration limits [as indicated in J29|, sec. II]] that 
specify the domain of feasible density matrices, directly within the Bloore-type framework, 
without transforming to partial correlations.) 

The computation of the m-th Hilbert-Schmidt moment over |p^^| is carried out in two 
stages. (The moments of the distribution of purity P{p) = Trp^ for quantum states with 
respect to tlie_ Bures measure have recently been determined by Osipov, Sommers and 



37 



Zyczkowski [36], while Giraud has investigated the Hilbert-Schmidt counterparts 

In the first stage, we perform an integration over the six-dimensional hypercube [—1, 1]^ of 

the m-th power of a (transformed) polynomial (P)-proportional to |/9^"^|-in seven variables 



29 



eq. (7)]). (The proportionality factor is (p22P33)^™-) The free (unintegrated) variable 
is of the form 



/^ 



PllP44 



P22P33 



where the p^j's are the diagonal entries of p. (In the related study [29|], 1/ = p^ was used as 
the principal variable, and in 25|], ^ = log p.) We have that (before the transformation ( TTOjl 
to partial correlations, yielding P) the multivariate polynomial (proportional to Ip^"'"]) 



P = -zf^fl'^ + 2Zu {Z12Z13 + Z2iZu) H^ + {V + W) I? + 2^23 (^12^24 + Zi'iZ^i) p 



^-23) 



(9) 



where 



and 



^ ~ (^34 ~" -'-) %2 ~ 2 (^14^23 + ^13^24) 2^34^12 + Z^^Z~ 



23 ^24 ^34 



w 



-2Z13Z14Z23Z24: + 2^13 (2^24 — 1) + 1. 



The transformation of the three correlations 213,2:14,224 to partial correlations (denoted 
2^13,2, 224,3, 2i4,23) takes the form [34] 



2;i4 -^ 212223234 + V2^2 - ly 2^23 - l2;i3,22;34 + Zi2\ Z^s - l\ zl^ - 1224,3+ (10) 



^12 ^ l\/^34 



'13,2 



1^/^24,3 



l2l4,23 + V ^12 ~ 1^23V 234 " 1^13, 22^24,3) 



2l3 -^ Z12Z2Z + \/ 2:12 ~ ly ^23 ~ 1^13,2) 224 ~^ 223234 + y 223 — 1^/2:34 — I224 



The jacobian for this transformation is (note that one of the six variables [or three partial 
correlations]-2;i4,23^s absent) 

-^^(2^12, 2^23, 234, 2i3,2, 224,3) = (212 - 1) (223 - 1) (234 - 1) Y 2:^3 2 - 1 Y2|4_3 - 1. (11) 



B. Intermediate functions/polynomials and their coefficients 

For the m-th moment {Momentm = (^), the indicated six-dimensional integration of 
pm ^^ j^Q^ reparameterized (partial correlation) form P^ over the hypercube defined by 
2^125 2:23, 2:34, -213,2, 224,3, 2i4,23 G ["1,1] takes the form-including a normalization factor of gf^- 
the ("intermediate function") result 

J{Zi2, Z23, ^34, 2i3,2, 2:24,3) [-^(2^12; ^235 2:34, 213,2, 2:24,3, 2i4,23)]™d2;i2d2;23d234d2;i3,2dZ24,3d2i4,23- 

For the first (m = 1) moment, we have the result (Fig. Hj) 

/,(„) = -^ + 2V-l. (13) 

^^' 5 125 5' ^ ^ 

for the second moment (m = 2), 

, , , 3/1^ 12/i6 20898/i^ 12/^2 3 
^^' 35 875 42875 875 35' ^ ' 

and for the third (m = 3), 

/ii2 54^10 27873/i^ 466876/1^ 27873/ 54/i2 1 
^^^'* ~ ~^ ~ "875 42875 1157625 ^ 42875 ~ 875 ~ 21' ^ ' 

At this point, we omit terms of lower order 2j than 2m, since their coefiicients-in a sym- 
metrical manner-match the coefficients C4m~2j{m)- Then, 

^ ^16 584/ii^ 278884^12 8984^1° 65788454/1^ 
^^^' ~ "33" ^ 5775 ^ 282975 ^ 4851 ^ 20543985 ^ ' ' ' ' ^ -^ 

3/x20_18/xi8 70881/xi*^ 2178728/xi^_ 59472398/1^2 ^4103383444/iio 
^^^''"^Tis" 143 49049 441441 4855851 273546273 ^"'' ^ ^ 

^,, /i24 2556/^22 5454/1^0 3359372/ii8 

/g(/i) = i— + ^— + i— + ^— + (18) 

^^^ 65 17875 2695 315315 ^ ^ 

3273117//^^ 597414184/1^^ 173821048732/1^2 
^ 86515 ^ 7872865 ^ 1771394625 ^ " ' ' 

10 



_^28 ^ 4298/i26 ^ 826637/i24 ^ 165865636yu22 ^ 71226035/i20^ 

^^^'* ~ ""85" " 27625 303875 8204625 722007 ^ '^ 
1947049760374yuis 93373201818911;ui6 33225665966177656^^^ 



6711055065 167776376625 48487372844625 
3/^32 6672/i30 12986136/i28 4250871568/i26 3319251741068;u24 

^^^'' ~ ^23" ^ 40375 ^ 3674125 ^ 121246125 ^ 14670781125 ^ ^ '' 
755365923834768^22 2024301386770232/i20 61510285844520752^1^ 3853435310162220966/1^6 

^ 826454003375 ^ 826454003375 ^ 14049718057375 "^ 724564031244625 ' ' ' 

and 

//36 9774^34 651051^32 8355664/i30 18384996780/^28 4848288282648/i26 

'^^^' ~ ^133 56525 145775 146965 39122083 1944597655 

(21) 
133915228926036yu24 61222919937476688/1^2 396008663496240078/^20 2103161056387491292/^1^ 

15026436425 2809943611475 10677785723605 47564681859695 

For the nine cases (m = 1, ..., 9) we have so far been able-with considerable computational 

expense-to explicitly compute, the coefficients of the corresponding 4m-degree even poly- 
nomials /m,(/i) are, as already indicated, symmetric-for reasons not immediately apparent- 
around the /i2™ term. 

1 . Formulas for the coefficients of the intermediate functions and their root and pole structure 

The constant terms (which equal the coefficients of the /i"^"^ term) of the intermediate 
functions, used in the computation of the moments of |p^"^|, are expressible as 

Additionally, the coefficients of the second and {Am — 2) terms are 

„, > „ IS 3(-l)"m(2.n(4m - 5) - 15) 

C,(m) = C,„,Mm) = ,oo(„_.)(„,^.)(„^|) - (23) 

Further, the coefficients of the fourth and {4m — 4) terms are 

3(-l)"m(2m(2m(2m(8m(6m - 7) + 155) - 13) - 1017) - 315) 



C4{m) = C4m~4{m) 



19600 (m - I) (m - i) (m + |) (m + |) 

(24) 



These results were obtained using the "rate", guessing program of C. Krattenthaler. Then 
further, M. Trott was able to obtain the result (using the FindSequenceFunction command 
of Mathematica, which searches for a possible rational form) 

CQ{m) = C4m-6{m) = (25) 

11 



(-l)™(m - l)m(4m(2m(2m(m(4m(20m(4m - 11) + 173) - 4303) + 4733) + 14911) - 9165) - 4725) 

529200 (m - |) (m - f ) (m - ^) (m + i) (m + |) ' 

From the formulas for these coefficients, it is clear that the numerator of the coefficient 
C2j{m){= C4m-2j(^)) of n"^^ is a polynomial of degree 3j, and the denominator is a poly- 
nomial of degree j + 2. (The denominators are very simple in structure ( l30|) -as evidenced 
above.). For j = 0, the roots are — | and — |, and as j increases by 1, an additional root 
1 larger in value than the previous smallest is added. Thus, poles occur at the coefficient 
functions at such half- integers.) 

2. Asymptotic convergence of dominant roots to half-integrs 

Utilizing this observation, we were then able to move on to obtaining the coefficients 
Cs{m) = C4m-s{m), Cio(m) = C4m~io{rn) Cuim) = Cim-uim), Cu{m) = CAra-u{m) 
and Ci6(m) = C4m-i6(^)^but not yet higher. In studying the (nontrivial) roots of these 
functions, we have detected one quite interesting feature. That is, as j increases, the dom- 
inant roots of C2j{m) show very strong evidence of converging to j — |, the subdominant 
roots to j — |,...For instance, for j = 8, the dominant roots of CiQ^m) = C4m-i6('T^) are 
7.49999796,6.4999352,5.4980028,4.4493216, while for j = 7, they are 6.5000204, 5.500556, 
4.515944. Such roots would then come increasingly close to canceling the near-to-matching 
poles in the denominators in C2j{'m) as j increases. (This behavior suggests an intimate 



connection with the theory of angular momentum and its [semiclassical] asymptotics 40|.) 



C. Use of the intermediate functions Imilj) to compute the m-th moment of \p^'^\ 

In the second stage of our procedure to compute the m-th Hilbert-Schmidt moment, we 
reverse the substitution (|8]) in these 4r?2-degree polynomials, multiply the result by the nec- 
essarily nonnegative factor (p22P33)^™ (the factor {p22Pz^Y had been removed, as previously 
indicated, in forming the polynomial P in seven variables, proportional to |p^"^|) and also 
by the jacobian corresponding to the transformation to Bloore (correlation) variables for the 
two-rebit density matrices p 

jac = {pnp22p33pM)^ , f3 = 3. (26) 
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{(3 = 6, for the two-qubit case, and 13 = 12 for the corresponding quaternionic scenario, in 
accordance with random matrix theory resuhs.) The result of the indicated multiphcations 
is, then, integrated over the unit three-dimensional simplex, 

Pll + P22 + P33 + P44 = 1, Pu>0, 2 = 1, . . . , 4 (27) 

to obtain the m-th moment (Cm)- I^ other words (taking into account the appropriate 
normalization factor), and setting P44 = 1 — pu — P22 — P33, the computation takes the form 

1146880 r^ r^-f^^ ri-Pii-P22 
Momentm = Cm = ^ / / / (28) 

TT^ Jo Jo Jo 



(P22P33)^'"(pllP22P33P44)2/m(W-^^^)dp33dp22dpll. 

V P22P33 

We are, in fact, able to perform the indicated symbolic integration, obtaining thereby (using 
an index i = 2j) 

Momentm = C = -^^^^^^t-0 2 4 ^ ('-^) T(-- + 2m+-) Q(m) (29) 
2293760 /_o^_o ^fi + 5V^f i „ 5^^ 



7r2r(4m + 10)V *='^'2'*- V 2 y V 2 



^lT.-o\T{^^] T{-^ + 2m+^\ Cdm)] + 

4 



1146880 ^ / 5 , ^ , , 

-r I m+ - I C2m{m), 



7r2r(4m + 10) V 2^ 
where the Ci{m){= C2j(wi))'s are our previously- indicated rational coefficient functions 
- (123]) ). symmetric about 2m. These (rational functions) Cj(m)'s themselves-as dis- 



cussed above-are the ratios of polynomials in m of degree ^ divided by the term (using 
the Pochhammer symbol, as well as rising factorials for gamma functions with half-integer 
arguments) 

,^, ,, fl-i \ r(m + |) 2^+2(2m + 3)!! , , 

denommator(a(m)) =[-^ + m) ^^^ = ^(,.+^+1) = (-, + 2r.-l)!! ^^'^ 

1-k^ 



nU-2,o....("^ 



2 ' 

For z = 4, by way of example, this gives us the denominator of fl24|) . 

On the other hand, the numerators of the Ci(m)'s for m > have zero as a trivial root, and 
for m > 4n, trivial roots 0, . . . n. 
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(2m + 3)!!)'c2™(m) + 2S2-o;|, (^(3 + 2)!!(3 - 2 + 4m)!!)'Q(m) 



Again, converting gamma functions with half-integer arguments to rising factorials, we 
have, equivalently to fl29|) . that 

Momentm = Cm = (32) 

r)7— 4m 

r[4m + 10] 

Pursuant to these formulas, the first moment (mean) of the Hilbert-Schmidt probability 
distribution of \p^'^\ over the interval [—jq, ^] was found to be 

Ci' = = ^ ^ -0.0011655, (33) 

^^ 858 2- 3- 11 -13 ^ ^ 

falling within the [negative] region (Ip'^^l < 0) of entanglement. (We were also able to 
compute this result using the alternative Euler-angle parameterization of the real density 
matrices 39| , but only after correcting a typographical error in the associated Haar measure 
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eq. (48)] , in which the factor sin X3 had to be replaced by its square. Also, we depart from 
the standard convention of denoting moments by /i, since that symbol has been employed in 
our earlier studies [25|, |29| and above (jH]).) Then, successively, the ([necessarily] decreasing 
in absolute value) raw (non-central) moments are 

^' = = ^ 0.0000108462, (34) 

^^ 2489344 2^0 ■ 11 ■ 13 ■ 17 ^ ^ 

C; = ^-^ = -^:^-, ^-^ ^ -1.2629773 ■ 10-^ (35) 

^^ 66216550400 2^3 . 52 ■ 7 ■ 11 ■ 13 ■ 17 ■ 19 ^ ^ 

C' = ^^^^Q = 21859 ^ 2 09311 • 10-^ (36) 
^^ 10443295948800 2^7 . 3 ■ 5^ ■ 11 • 13 • 17 ■ 19 • 23 ' ' ^ ^ 

23071 23071 „ , , 

G = = — r. -0 ~ -4.27531 • 10^^\ 37 

^^ 539633583390720 2i8 • 3 • 5 ■ 7^ ■ 13 • 17 ■ 19 ■ 23 • 29 ' v ; 

C' = ^^ll^^l = 7-43-103-107 ^ ^ Q^g^g _ ^^_,, 

^'^ 3253917653076541440 228 ■ 3 ■ 5 ■ II2 ■ 17 ■ 19 ■ 23 ■ 29 ■ 31 ' ' ^ ^ 

C' = 419856257 ^ 

^^ 15366774022001834065920 ^ ^ 



43 • 2179 ■ 4481 



230 . 34 ■ 5 ■ 11 ■ 13 ■ 17 ■ 19 ■ 23 ■ 29 ■ 31 ■ 37 



-2.73223 ■ 10 



-14 



16945249 
^^ ~ 21117403549591928832000 ~ ^ ■^ 



109 ■ 155461 



8.02431 ■ 10 



-16 



233 ■3- 53- 11- 19. 23 -29 -31 -37 -41 
and (requiring four days of Mathematica computation on a MacMini machine) 

,_ 6102620963 

^^ ~ "240565904621616585139814400 ~ ^ ^ 
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19 ■ 199 ■ 1614023 



-2.53678 ■ 10 
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237 . 3 ■ 52 ■ 113 . 13 ■ 23 • 29 • 31 ■ 37 ■ 41 ■ 43 
(After four weeks of uninterrupted computation, we did not succeed, however, in determining 

CioO 

Interestingly, the sequence of denominators immediately above (in apparent contrast to 
that of the numerators) appears to be "nice" in that the number of their prime factors do 
not grow rapidly, but rather linearly. This is a strong indication of the possible existence of 



a "closed 
factorials 



brm" , that is an expression which is built by forming products and quotients of 



41 



fn. 12]. 



D. Use of moments to estimate the probability distribution over \p 



PT\ 



In Fig. O we display a (naive) fit of a simple power series in \p\^'^ of degree nine to the 
computed first nine moments ( fl5^ - fHT]) ) (cf. [37|, Figs. 1, 2]) of the Hilbert-Schmidt proba- 
bility distribution over |p^"^|, where p is a generic two-rebit density matrix. No nonnegativity 
constraints were, however, imposed and considerable incursions into negative regions result. 
(Such negativity can be obviated through the use of maximum-entropy, spline-fitting and 



other methodologies [27, 



28 



30 



m 



46 



471], and we do explore such directions.) 



The Hilbert-Schmidt separability probability predicted by the curve in Fig. |5]-that is the 
"probability mass" (the resultant of both positive and negative values) lying within the 



interval [0, 2^6]^^^ 0.39648, while our previous studies [2^], indicate that the actual value is 
somewhat higher, k. 0.45-a discrepancy the use of additional higher-order moments should 
ameliorate. 

Since the plotted distribution (Figj5]) appears to be unimodal, one can presumably use the 
computations of the first and second moments above to isolate the mode of the distribution 



within the interval 



43 



eq. (13)] 



30397 

51 



30397 



51 



858 



4576 4576 



858 



{-0.00650062, 0.00416962}, 



(42) 



containing the value |p^^| = 0. Narrower intervals containing the mode 
using higher-order moments and the associated Hankel determinants 43 



can be obtained 
Thm. 3.2]. 
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E. Numerical computations of higher-order intermediate functions (j > 9) 

It is clear that it would be of considerable utility to have available exact values for still 
higher-order (than m = 9) moments-and for the coefficients C2j{m) = Ci{m) of the terms in 
the intermediate functions/polynomials Im{fJ')-hut the associated computational demands 
seem quite considerable. 

Our only current recourse, in this regard, then, appeared to be a numerical one. 



We employed a quasi- Monte Carlo (Tezuka-Faure [48|, |49|) procedure (using 18,870,000 
[low-discrepancy] points) to estimate the values of C2j{m) for j > 9 and m > 10, for 
m = 10, ... 50. (Actually, by the evident symmetry around the 2m-th power of the coeffi- 
cients of the intermediate functions, we were able to obtain two values [which we averaged] 
for each point.) Coupling these approximate results with the exact formulas for the inter- 



mediate functions obtained above for j < 9 (sec. IIIB|) . we obtained estimates of the first 



fifty moments. (We investigated the possibility of using these additional numerical results 
to infer the desired further exact formulas for the C2j(m)'s, but this seemed too demanding 
a task, given the ordinary machine precision of our quasi- Monte Carlo simulations and the 
evident high [multi-digit] complexity of the exact coefficients.) 

1. Cumulative distribution function calculations 

We then used these values (the first nine moments being exact, and the remaining "semi- 
exact" forty-one, being sums of exact and numerical terms) in the procedure of Mnatsakanov 
for approximating the "moment-determinate cumulative distribution function (cdf) from its 
moments". The relevant formula for the cdf (at the separable/nonseparable boundary of 
principal interest to us) based on the first K = m moments (linearly mapping \p^'^\ to lie 



in [0,1], and the moments accordingly) is of the form |27i, eq. (2)] [50|, eq. (1.3)] 

(We have to subtract this [entangled probability estimate] from 1 to get the [complementary] 
separability probability estimate, to be plotted as a function of the number of moments K.) 
In Fig. [6], we show estimates of this value based on increasingly large numbers of moments. 
(For more than thirty-six moments, the results turn negative. If we just employ in this same 
reconstruction procedure, the known exact formulas for the coefficients of the intermediate 
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functions-effectively setting the supplementary/corrective numerical terms to zero-then the 
two sets of results are quite close up to twenty moments, but then become considerably more 
ill-behaved when none of our supplementary numerical results is included [Fig. [7]. In a recent 
interesting study, Gzyl and Tagliani concluded that thirty-two was "the maximum allowable 
moments before incurring numerical instability, unless one conducts the calculation with 



high accuracy [50|]. They also assert that "the additional information introduced by using 
the (M + l)-order moment is Visible' only after the 0.6M-th decimal digit". In this regard, 
let us note that the quasi-Monte Carlo calculations used to complement the exact results 
here were conducted with only ordinary machine precision.) 

It would appear that the wide range and lack of stability of estimates is reflective of the 
ill-posedness of the Hausdorff moment (inverse) problem (stemming ultimately from the lack 
of orthogonality of the sequence 1, x, . . . , x", . . . [50|). In Fig. |8] we show the ratio of the 
exact (but incomplete for m > 9) moment computations to that (semi-exact one) based on 
the exact and complementary numerical results. 

F. Libby-Novick (three-parameter) probability distribution 



We have explored the use of extensions of the (two-parameter) beta distribution 5l|, chap. 



5] 52| to better fit the nine exact moments than found, as presented above (sec. IIP 1\\ . with 
the particular distribution ([5]) that had been fit to the first two moments of |p^"^|. Doing so, 
we were able to obtain a t/iree-parameter Libby-Novick (LN) distribution of the form [53] 
d, eq. (IX. 1)], 

a ^ 3.7141606, b ^ 359.577737, A ^ 0.00064805. 

This gave us a further considerably improved fit over that of (jS]) to the first nine exact mo- 
ments. (The ninth moment was now predicted within 99.6%-as opposed to within 98.6%-and 
the preceding moments better still-in a monotonically declining goodness-of-fit manner from 
the first to the ninth.) The estimated separability probability now increased to 0.429121. 

To arrive at the probabihty distribution ( l44l) . we started with the beta distribution fit ([5]) 
that exactly reproduced the first two moments, now trying to fit the first three moments. This 
involved a very long (slowly converging) iterative process, which at each stage, appealingly. 
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seemed to improve the fit to all our computed fifty (nine exact and forty-one, "semi-exact") 
moments. Additionally, the estimated separability probability seemed to increase at each 
step, which we found to better accord with our earlier extensive numerical investigations 
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, sec. IX. A] |29|, sec. V.A.2] [25]. 



One can employ the Libby-Novick distribution flH|) as a "baseline density", in the manner 



28 



eq. 



below using the beta distribution, following the methodology of Ha and Provost 
(3.5)], to generate estimates of the HS separability probability. 

G. Moments based on lesser principal minors of p than its determinant 

A necessary, but not sufficient condition that a two-qubit density matrix be separable is 



that any 3x3 principal minor of its partial transpose be nonnegative 25| . So, we can select 
one such minor, say (cf. (jH])), 

2 

minorsxs = "/'^ (/x^2:i4 - 2^Zi2ZviZii + zl^ + zl^ - 1), (45) 

/i 

expressed in terms of the Bloore parameterization, and compute the associated moments. 

(Enforcing the nonnegativity of such a minor yields an upper bound of ||_~ 0.628571 on the 



Hilbert-Schmidt probability of separability of generic two-rebit systems 25|.) We have been 



able to compute exactly the first forty-five such moments of the probability distribution 
over the interval [— |i ^]5 following very much the same scheme as we pursued for the first 
nine moments of \p^'^\- The first two moments are —i^ and -j^^, while remarkably, the 
(raw) third moment is identically zero (cf. [54j). The form the corresponding "intermediate 
functions" now took for these first three cases were (cf. flT3|) -f fT5|) ) 

/i(^) = ^(^'-3), (46) 

/2(;u) = ^(75/x^-182//2 + 395) 

and 

^ , , 125^^ - 297^^ + 675/i2 - 935 

^'^^^ = 2625 • 

General formulas could now also be obtained for the coefficients of the 2j-th powers of 

the intermediate functions. These involved the Mathematica functions LerchPhi or Differ- 

enceRoot (cf. (!22|) - (l25|) ). However, we were unable to find a general formula for the m-th 

moment, even based on the availability of the first fifty-eight moments. 
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III. MOMENTS OF THE PRODUCT OF |p|^^ AND \p\ 

The main foci of our initial analyses had been, firstly, the moments of the determinant 
Ip'^^I of the partial transpose of generic two-rebit density matrices, and, secondarily, the 
moments of the (necessarily nonnegative) determinant |p| of the underlying density matrix 
itself-the moments of the latter determinant being more computationally amenable, it is 
clear, to exact analyses. 



Now, in this context, it is not unnatural to ask the question of t 
of the product of these two determinants, |p||p^"^| = \pp^'^\ (cf. 55 



le nature of the moments 
p. 564]). This approach 



potentially contributes insight into the separability probability question, since a value of 

the product less than zero still indicates the presence of an entangled state, and a value 

n 

greater than zero, a separable state. (It would appear that the argument in [3] of Augusiak, 
Horodecki and Demianowicz could be adapted, so that the value of the product could be 
construed as that obtainable from a single certain observable measurement.) 

In undertaking the associated analysis of moments, we immediately encountered a most 
interesting result. The first moment or mean of the normalized product |p||p^"^| is zero, that 
is 

C; = 0, (47) 

the associated "intermediate function" being (cf. ( IT3l) ) 



(Numerically, the first moment of the absolute value of the normalized product is 5.86519 • 
10~^. This is close to the absolute value, i^^^i^^ ~ 5.09397- 10^^, of the product of the mean 
of |p|, that is, 2^ and the mean of |p^"^|, —gig-) The second moment (cf. (IMll ) is 

Co = = -r. — -. — 7, ^ 1.2288585 - IQ-^^ ^49) 

^ 5696343244800 2^^ -3^-5^-13-17-19-23 ^ ' 



(so, the corresponding standard deviation is the square root of this, that is -j^^^l'' ~ 1.10854- 
10~^) with the associated intermediate function being (cf. f lT^ ) 

192/ _ 12032//6 5561984/i^ _ 12032^2 192 
'^^^' ~ 94325 " 1528065 ^ 184895865 ~ 1528065 ^ 94325' ^ ^ 

(To compute the m-\ky moment of the probability distribution of the product, using the new 
set of intermediate functions, we employ the same formula as (|28|) . but for the replacement 
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of the exponent | by | + m.) We see that the coefficients of the constant (and highest power 
of fi) terms in the first and second intermediate functions immediately above (dHj), fl50|) ) are 
— ^ and gillg • The next four such constant coefficients (m = 3, 4, 5, 6) have been found to 
be -iiis' 586831245 ' " 967™577 ^nd 163^11875 - However, we have yet to obtain a general 
formula, parallel to (1221) . in light of the increased computational burden, encompassing these 
six values. 

The third to the sixth moments are 

c;- 



677899511057612800 2^^ • 32 • 52 • 13 • 17 • 19 • 23 
1 7 



45973294808920227840000 21^ • 32 • 52 • 13 • 17 • 19 • 23 

1 7 



1.2288585 • 10"^^ 


(51) 


^ 1.2288585 ■ 10"^^ 


(52) 



1-12 



/• __ __ r~^ 1 9288^85 10" 

^^ 11662680803407302839532257280 2i8 ■ 32 ■ 52 ■ 13 ■ 17 • 19 • 23 '^ 

(53) 

and 

3929 7 

Cfi = = -r. — -0 — -0 ~ 1.2288585-10 

^'^ 4158654163938276392103553381781471232 2^» ■ 32 ■ 52 ■ 13 ■ 17 ■ 19 • 23 

(54) 
Despite their lengthy digital descriptions, the ratios of these six moments to the HS 
moments of Ipp'^-given by ([2]) are rather remarkably simple, that is, 

77 24 209 598 3929 
^' 54' 55' 175' 833' 3724^ ^^ 

We have only so far been able to compute the two-qubit analogue of the first (zero) of the 
six ratios above-it turning out quite remarkably to be |/ Since these ratios are so simple, it 
suggested to us that we might be more able to progress in our series of analyses, by making 
our initial goal the computation of these (unknown) but apparently well-behaved ratios-for 
higher-order moments-rather than the very small values of these moments themselves. 

IV. SUMMARY 

We have studied here the moments of probability distributions generated by certain de- 
terminantal functions of generic two-qubit density matrices (p) with real entries ("rebits") 
over the associated nine- dimensional convex domain, assigned Hilbert-Schmidt measure. It 
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was found that the mean of the (nonnegative) determinant \p\ is 2^, the mean of the deter- 
minant of the partial transpose Ip'^"^ [-negative values indicating entanglement-is — g|g, while 
the mean of the product of these two determinants, IpHp"^^! = |pp^"^|, is zero. We deter- 
mined the exact values-also rational numbers-of the succeeding eight moments of Ip'^"^]. At 
intermediate steps in the derivation of the m-th moment of Ip'^'^l, rational functions C2j{'m) 
emerge, yielding the coefficients of the 2j-th. power of even polynomials ("intermediate func- 
tions" Imi/J')) of total degree 4m. These functions possess poles at finite series of consecutive 
half-integers {m = — |, — |, . . . , ^^f^), and certain (trivial) roots at finite series of consecutive 
natural numbers (m = 0, . . . , |_^J)- The (nontrivial) dominant roots of C2j{jn) appear to 
converge, as j increases, to the same half-integer values {m = . . . , -^^, "^O- If formulas 
for C2J {tti) could be developed for arbitrary j-we do possess them already for j < 9-then, 
the desired Hilbert-Schmidt separability probability would be computable to high accuracy. 
We reproduced the (linearly transformed) first nine moments of |p^^| quite closely by 
a certain (two-parameter) beta distribution, and still more closely by a three-parameter 
(Libby-Novick) extension of it. The first two moments of |p^"^|-when employed in the one- 
sided Chebyshev inequality-gave an upper bound of |^||| ~ 0.874759 on the Hilbert-Schmidt 
separability probability of two-rebit density matrices. We ascertained by numerical methods 
that the orthogonality established of \p\ and \p^'^\ with respect to Hilbert-Schmidt measure 
does not hold with respect to the Bures (minimal monotone) measure, nor if we slightly 
distort the symmetry of our basic nine-dimensional generic two-rebit scenario. 
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FIG. 2: Ratios of the m-th. exact moment of \p \ (linearly transformed to the interval [0,1]) to 
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FIG. 6: Separability probability estimates based on differing numbers (exact for r?i < 9 and "semi- 
exact" for ?n, > 9) of moments, using the reconstruction procedure of Mnatsakanov [27|. The 
horizontal a xis is drawn to intercept the vertical at ^^ w 0.537619, the least upper bound so far 
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FIG. 7: Separability probability estimates, without any complementary numerical (quasi- 
MonteCarlo) input (for m > 9), using the reconstruction procedure of Mnatsakanov [271]. The 
^lorizontal axis intercepts the vertical at ^^ ss 0.537619, the least upper bound so far established 
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FIG. 8: Ratios of moments computed using the exact and complementary numerical (for m > 9) 
results to those based only on the exact, but incomplete (for ?7i > 9) results 
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